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ABSTRACT 



We examine the observable properties of simulated barred galaxies including ra- 
dial mass profiles, edge-on structure and kinematics, bar lengths and pattern speed 
' evolution for detailed comparison to real systems. We have run several simulations in 

which bars are created through inherent instabilities in self-consistent simulations of 
a realistic disc+halo galaxy model with a disc-dominated, flat rotation curve. These 
simulations were run at high (N=20M particles) and low (N=500K) resolution to 
0^ , test numerical convergence. We determine the pattern speeds in simulations directly 

from the phase angle of the bar versus time and the Tremaine- Weinberg method. 
Fundamental dynamics do not change between the high and low resolution, suggest- 
. ing that convergence has been reached in this case. We find the higher resolution is 

' needed to simulate structural and kinematic properties accurately. The edge-on view 

| of the higher-resolution system clearly shows the bending instability and formation of 

a peanut-shaped bulge. We determined bar lengths by different means to determine 
i-C ' the simulated bar is fast, with a corotation to bar length ratio under 1.5. Simulated 

bars in these models form with pattern speeds slower than those observed and slow 
q ■ down during their evolution. Dynamical friction between the bar and dark halo is 

5-h ' responsible for this deceleration, as revealed by the transfer of angular momentum 

c/3 , between the disc and the halo. However, even though the pattern speed is reduced at 

later times, the instantaneous scale length of the disc has grown sufficiently for the bar 
motion to agree with many observations. By using a different model and simulation 
technique than other authors, we are able to compare the robustness of these meth- 
ods. An ani mation of the face-on and edge-on view s of the 20M particle simulation is 



available at http : / / www, astro . ut oronto . ca^oncill 
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1 INTRODUCTION 

The dichotomy between the regular and the barred spirals 
in Hubble's original galaxy classification scheme is a long 
standing problem in galaxy evolution and dynamics. It is 
unclear what causes a spiral galaxy to become barred or not, 
or more quantitatively, why barred galaxies represent about 
30 to 70 per cent (deVaucouleurs 1963; including weak bars, 
Sellwood & Wilkinson 1993; viewed in IR, Eskridge et al. 
2000) of nearby galaxies. This picture becomes even more 
confusing at high redshifts, where Abraham et al. (1999) 
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have found the fraction to drop beyond a redshift of 0.5. 
Bars undoubtedly form from dynamical instabilities inher- 
ent in self-gravitating axisymmetric discs. Early simulations 
(e.g. Ostriker & Peebles 1973) have shown that purely self- 
gravitating discs are instantly susceptible to a bar insta- 
bility, while the addition of a surrounding spheroidal grav- 
itational potential can prevent bar formation; a more re- 
cent study suggests this is merely a slow down of formation, 
and the bar is actually enhanced at later stages (Athanas- 
soula 2002a). The surrounding stabilizing spheroids can be 
identified with the dark haloes believed to surround spirals 
and are responsible for the nearly flat rotation curves of all 
galaxies. Cold dark matter models of cosmology predict that 
most spirals are embedded within nearly isothermal haloes 
(Dubinski & Carlberg 1991; Navarro, Frenk & White 1996, 
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1997). It is therefore useful to examine the consequences of 
dark haloes on the formation and evolution of barred galax- 
ies to seek consistency with the prevailing world model and 
the observable universe. 

The pattern speed of the bar is an important although 
difficult to measure indicator of disc dynamics and dark 
halo structure. Observations of early-type barred galaxies 
(Merrifield & Kuijken 1995 (M&K); Gerssen et al. 1999) 
through the application of the Tremaine- Weinberg method 
(Tremaine & Weinberg 1984b) (T&W) have shown the pat- 
tern speed to be 'fast', since the bars end near corotation. N- 
body models of disc dynamics with static background haloes 
(summarized in Sellwood 1981) have shown the pattern 
speed also to be fast. However, soon afterward it was rec- 
ognized (Tremaine & Weinberg 1984a; Weinberg 1985) that 
rotating bars would spin down significantly due to dynami- 
cal friction from the dark halo. More recent simulations with 
self-consistent disc and halo distributions clearly showed this 
effect: angular momentum is transferred from the bar to the 
halo through torques due to gravitational wakes in the halo, 
which result in surprisingly low pattern speeds for simu- 
lated bars (Debattista & Sellwood 2000 (D&S); Misiriotis 
& Athanassoula 2000). D&S conducted parametric studies 
by varying the halo to disc mass ratio of their models and 
found that a maximal disc yields the least pattern speed 
slow down. They argued that the observational evidence of 
fast bars implies a maximum disc in all barred galaxies. 

Bar structure and dynamics can also be examined in 
edge-on systems. The thin bars that form in discs are sub- 
ject to a buckling instability which causes a bar to bend 
vertically and thicken into a bulge-like object within a few 
dynamical times, as was shown numerically by Raha et al. 
(1991), Pfenniger & Friedli (1991), and Combes & Sanders 
(1981). Kuijken & Merrifield (1995) (K&M) have shown that 
the kinematics of peanut-shaped bulges of NGC 5746 and 
NGC 5965 are consistent with orbits in a bar potential, the 
theory of which was confirmed by Bureau & Athanassoula 
(1999) (B&A). The observational test described by K&M 
can also be applied to simulations, which can be viewed at 
any angle for more complete results. 

Through B and I band observations of 15 galaxies, 
Elmcgrccn & Elmegreen (1985) have shown that early- (SB0, 
SBa) and late-type bars have significantly different proper- 
ties. Although resonance positions and bar lengths are hard 
to quantify observationally, there is evidence to suggest that 
early-type bars are longer, out to the corotation radius, have 
a flat surface brightness profile along their length, and end 
beyond the turnover radius of the rotation curve. The late- 
type bars on the other hand are much shorter, with expo- 
nential surface brightness profiles. 

In this paper, we revisit the problem of the bar and 
buckling instabilities in self-consistent disc galaxy models 
with live haloes. Our goal is to study the pattern speed evo- 
lution of numerical bars as well as quantify the structural 
and kinematic properties for comparison with observed face- 
on and edge-on barred galaxies. We have taken care to re- 
scale our models to observed systems for comparison. We 
extend previous work by introducing a new, more realis- 
tic mass model and going to much higher resolution with a 
simulation containing N=20M particles, a factor of 20 times 
larger than most work on the subject. Our motivation for 
using high resolution is to examine the numerical conver- 



gence of results, since subtleties of dynamical interactions 
between discs and haloes may not be captured correctly even 
with haloes with more than « 1M particles (e.g. Weinberg 
2001). The disc of our large simulation contains 10M parti- 
cles, taking us well out of the regime where disc self-heating 
contaminates results. We are confident that the dynamics of 
these models reliably represents the gravitational behaviour 
of real galaxies. There are also disagreements between var- 
ious simulations of bars in the literature which use differ- 
ent N-body codes and initial conditions. We attempt here 
to achieve convergence by comparing three of our simula- 
tions with those of other groups, as well as to observations, 
by examining the bar structure, pattern speed, and edge-on 
kinematics. 



2 MASS MODELS 

We simulate bars by setting up an initially axisymmetric 
system including a disc and dark halo which is formally in 
equilibrium but unstable. Models with rotation curves which 
are disc-dominated, approaching the maximal disc as defined 
by van Albada & Sancisi (1986), usually form a bar within 
a few dynamical times, and the Ostriker & Peebles (1973) 
empirical criterion is still a useful indicator of the instability. 

We use the method of Kuijken & Dubinski (KD) (1995) 
to generate self-consistent disc+halo models with a nearly 
flat rotation curve. We consider models with disc-dominated 
rotation curves with and without a central bulge. The KD 
models are generated from a distribution function (DF) that 
is the sum of up to 3 functions: a three-integral disk DF, a 
bulge DF modelled as King model with an energy cut-off 
E < and a halo DF that is a flattened King model DF (or 
lowered Evans model) with the usual truncation at a tidal 
radius at E — 0. We examine models with and without a 
bulge as described below. The disk DF is a function of E, 
z-angular momentum, L z and a third approximate integral, 
the z energy, E z = i 2 /2 + &(z). The disk DF is constructed 
assuming a exponential radial surface density profile and 
an approximately sech 2 z edge-on profile with fixed vertical 
scale length Zd- The disk squared radial velocity dispersion, 
crfj is assumed to decline exponentially with the same scale- 
length as the disk like real galaxies. We can generate N- 
body realizations of these models and they are formally in 
equilibrium with an initial virial ratio 2T/W = —1.0. 

The bulgeless model is nearly a formal maximal disc 
model, with the rotation curve rising to a roughly flat pro- 
file within two scale lengths (Figure 1). Although the halo 
DF is a King model, the mass profile does not have the usual 
core since the disk mass dominates the centre. The halo den- 
sity profile has a mild cusp with p cx r~ (,,T to within 0.1 
disk exponential scale lengths similar to an NFW profile. 
Whether or not the haloes of spiral galaxies have a central 
core is still controversial but our model halo here is consis- 
tent with rotation curve decompositions (e.g. Kent 1986). 
The disc is in formal equilibrium with a Toomre Q ~1.1 
measured at R=2.QR A and is approximately constant in the 
range 1.0 < R/Rd < 4.0 rising to higher values beyond 
these limits. While the model is stable against axisymmet- 
ric perturbations it is inherently unstable to bar formation 
because of the dominance of the disk in the rotation curve. 
The model is scaled such that G=Rd=v max = l , and the to- 
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Figure 1. Rotation curves for disc+halo (left) and disc+bulgc+halo (right) models in simulation units. The contribution of the disc is 
shown in dots, the halo in long dashes, the bulge in short dashes and the total as a solid line. R is measured in exponential disk scale 
lengths. 
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Figure 2. Initial halo density profile: a King model halo with a 
very small core; the profile is similar to an NFW profile with a 
-0.7 power law slope until at least R ~0.1 

tal disc+halo mass is about 12 (disc mass of 1.98 and halo 
mass of 10.57). We realize these models as N-body systems 
at high resolution with N=20M particles (10M each in the 
disc and the halo) and low resolution with N=500K parti- 
cles (250K each in the disc and halo). Thus, all quoted sim- 
ulation lengths are in terms of the exponential disc length, 
and all velocities are fractions of the maximum. We also 
consider a N=500K model with a small compact bulge of 
mass M — 0.32 with halo and disk mass profiles essen- 
tially unchanged. The Toomre Q also remains comparable at 
~ 1.06. The code and parameters for generating these mod- 
els are available by request. Our models differ from D&S who 
use a Kuz'min- Toomre disc with a vertical Gaussian profile 
and a lowered polytrope distribution function (DF) for the 
halo. The recent publication by Athanassoula & Misiriotis 
(2002) (A&M) contains models based on Hernquist's (1993) 
method. 

As well as specifying the following results in the origi- 
nal simulation units, the appropriate lengths and times are 
scaled to match the SBa galaxy NGC 4596 for easier com- 
parison with real galaxies. We use the scale length of 3.2 kpc 
and rotational velocity of 120 km s _1 given by Gerssen et 
al. (1999). The sech 2 scale height for the disc is 0.1 simula- 



tion units or 320 pc. For these scales the unit of time is 26.3 
Myr. For the runs, we use a time step of At = 0.05 (1.32 
Myr) for 10,000 time steps, for a total t=500 (13.2 Gyr) 
or approximately the Hubble time. We soften gravitational 
forces using a Plummer law with softening length e = 0.0125 
(40 pc) for the halo and e = 0.005 (16 pc) for the disc. All 
runs are executed with a parallel tree code (Dubinski 1996) 
running on a PC cluster or a 32 processor Compaq GS320. 
The 20M particles simulations needed 2 minutes per step on 
the GS320. The total energy drifts by no more than 1 per 
cent, and total angular momentum is conserved to typically 
within 0.1 per cent. The results discussed below are for the 
high resolution simulation unless otherwise noted. 



3 BAR FORMATION AND STRUCTURAL 
EVOLUTION 

Although these models are formally in equilibrium they are 
strongly unstable. A bar develops by t=26 (686 Myr), and 
is sustained for the duration of the run. Visually, the bar 
grows to its greatest extent at around t=50 (1.32 Gyr). The 
bar then bends and buckles between i=70 and i=90 (a pe- 
riod of 526 Myr) to leave a slightly slower bar of comparable 
length. There are isophotal twists seen during the bar buck- 
ling phase, but as noted in Shaw et al. (1993) these disappear 
with the increased heating of the disc seen during the bar 
evolution. After this short phase, the bar shows elliptical 
isophotes with a butterfly- or dumbbell-shape when viewed 
face-on. Fig. 3 shows the face-on views for t=50 to 500. 

The evolution of the bar viewed edge-on clearly reveals 
the buckling instability. Initially the disc is thin; the for- 
mation and motion of the bar heat the disc vertically and 
azimuthally. At the onset of buckling there is a very short 
lived phase of about 260 Myr in which the bar bends into 
an arc before buckling into a boxy shape object (edge-on 
views are shown in Figs 4 and 5). The boxy-shaped bulge 
remains after buckling, and slowly evolves into the famil- 
iar peanut shape, which is most prominent when the bar is 
viewed edge-on, with its long side in the plane of the sky. 
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appigx. 1 scale length 

Figure 4. Logscale density diagrams of the edge-on views of the simulated galaxy at t=50, 76, 100, 250, 350, and 500. 




' scale length 

Figure 5. The edge-on views of the simulated galaxy showing the warp at t=80, and the end on views of the bar at t=150 and 500. 
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Figure 6. Rotation curve of the simulated galaxy at t=0, 100, 
200, and 400. The bar causes an increase in central mass, which 
increases the initial slope of the rotation curve as well as raising 
the maximum value. 



When viewed end-on, the galaxy resembles a disc with a 
spherical bulge. The bar remains an obvious component of 
the disc and shows no sign of collapsing into a spherical 
bulge. There is an animation of both face-on and edge-on 
views available at www.astro.utoronto.ca/ oneill. 

There is an increase in central density following the ini- 
tiation of the bar, which influences the initial slope of the 
rotation curve. Fig. 6 shows the rotation curve at t~0, 100, 
200, and 400, and it shows the increased slope and peak dur- 
ing the simulation. The surface density along the bar is flat 
at small R for t=50, but rapidly turns exponential following 
buckling. These data will be discussed in more detail in the 
following sections. 

Even without the gaseous component, this simulated 
galaxy can be compared with observations, and some obvi- 
ous features stand out. The face-on views look very much like 
observed bars, and the edge-on views are similar to several 
studies comparing peanut-shaped bulges and bars (K&M; 
Merrifield & Kuijken 1999). More detailed comparison with 
the latter follows in Section 5. 

Recently, A&M conducted a series of tests to compare 
N-body simulations of different mass models, which are ide- 
ally suited to comparison with our results here. They sim- 
ulated three different models: one disc-dominated (MD), 
which is similar to ours here, another disc-dominated model 
with an added bulge (MDB), and one halo-dominated (MH): 
all created using the Hernquist (1993) models. 

Our work agrees quite well with A&M on isophote shape 
and rotation curve progression for the first half of the simu- 
lations, which span approximately the same time. Our evo- 
lution here disagrees to some extent with their disc domi- 
nated model for the latter half of the run: the mass build 
up at the centre of the disc causes the inner slope of the ro- 
tation curve to become more steep (albeit slightly) and the 
maximum velocity to increase even during the final stages; 
whereas, A&M found little difference between the rotation 



curves for the last quarter of their run. This also carries 
over to the bar isophote structure, as our increasing cen- 
tral density causes the bar shape to become similar to their 
MDB model. Therefore, the latter half of our simulation 
is more readily compared to dumbbell-shaped, face-on and 
peanut-shaped, edge-on isophotes, which they produce with 
the MDB model. Our low-resolution model, with a similar 
number of disc particles as were used by A&M, shows only 
boxy isophotes, but was run for half as long as the high- 
resolution model. Since the buckling stage was later in this 
simulation, the evolution into a peanut-shaped bulge would 
be expected after the end of the simulation, so it is uncer- 
tain whether we would agree with A&M in this case or not. 
The rotation curve for our low-resolution model does not 
change slope for the last 100 simulation time units, which 
may suggest that the mass redistribution does not continue 
as long as the high-resolution model run. 

We also investigated the evolution of the halo. Unlike 
Weinberg & Katz (2002), the central density distribution of 
our halo remained virtually unchanged. There were several 
differences between our models which may account for this. 
(1) Their bars were much stronger than ours (30 per cent 
of the disk mass), appearing earlier in the formation of the 
galaxy out of a cold, gaseous disc (ours formed in the stellar 
disc of a galaxy initially in equilibrium), which will have a 
greater dynamical friction effect on the halo. (2) Our halo 
is based on the King model, which according to Weinberg 
& Katz (2002) is missing the key resonance for redistribut- 
ing the halo mass. However, the KD models are composite 
DFs and the resulting halo has a mild cusp with an effective 
power law slope of -0.7 at the centre, which is similar to an 
NFW profile. (3) They stress the need for high resolution, 
using 4M particles in the halo; they are concerned that tree 
codes induce more small scale noise and require higher res- 
olution - they did not outline how much more, but the 10M 
particles in our halo should be enough to see some indication 
of halo density redistribution. 



4 BAR PATTERN SPEED 

T&W describe the only method available to measure pattern 
speeds of real galaxies from observable quantities alone. This 
method uses luminosity, line-of-sight velocity, and tilt angle 
to determine the rotation rate of the bar. To date there 
have been a few early-type galaxy bars measured this way, 
and every one of them has been found to be fast, with the 
corotation distance to bar length ratio between 1 and 1.5. 
This result has been difficult to reproduce in simulations, 
as all the bars are greatly slowed due to dynamical friction 
between the bar and the halo. The amount of mass in the 
halo compared with the disc is still a point of contention, 
because the existence of a flat rotation curve constrains, but 
does not pinpoint, this ratio. Studies on this relation refer 
to 'maximal' or 'sub-maximal' discs, where the former has 
the largest disc to halo mass ratio in the inner regions of the 
disc that is consistent with observed rotation curves (van 
Albada & Sancisi 1986). D&S used the dynamical friction 
between the bar and halo to show that only a maximal disc 
is likely in barred galaxies, since these bars slow down the 
least. The final state of their maximal disc simulation was 
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Figure 7. Pattern speed of the bar vs. simulation time units 
(rad t -1 ): the pattern speed continues to slow even after a Hubble 
time. 

nearly a fast bar, with a corotation to bar length ratio of 
1.6 ±0.3. 

We used two different methods to determine the pattern 
speed of the bar (O p ). The first was a direct measurement 
of bar rotation by determining the orientation angle of the 
bar as a function of time. Fig. 7 is a plot of pattern speed 
versus time in the original simulation units. Early on, the 
pattern speed is 0.44 at t=26 but there is a rapid reduction 
during the initial bar growth, and by £=50 it has slowed to 
0.35, then there is another drop during the buckling phase 
from t=65 to 95, with the pattern speed settling to 0.33 (the 
following section compares these with observations). After 
that initial phase of rapid evolution, the pattern speed is in- 
fluenced by dynamical friction, as it gradually slows to less 
than 0.21 by the end of the simulation. Even after approxi- 
mately a Hubble time, it does not reach a steady state; it is 
still transferring angular momentum to the halo (see Fig. 8). 
We note in Fig. 7 that the pattern speed shows significant 
scatter early on, during bar formation and buckling phases. 
After formation, the pattern speed smoothly, but quickly, 
decreases. During buckling, the scatter is probably due to 
difficulty in precisely measuring bar position in the midst 
of twisted isophotes. There is also a series of oscillations 
which continue after the buckling phase which we believe to 
be real, whose period corresponds to slightly more than the 
circular period of the disc at corotation. We are uncertain 
at this time of the cause of these oscillations. 

In order to compare these results with both observations 
and other simulations, they must be scaled to comparable 
units. We have taken the R (scale length) divided by v (cir- 
cular velocity) to determine a scaling factor (DS subscripts 
for Debattista & Sellwood, OD subscripts for O'Neill & Du- 
binski): 

QpDS _ vds x Rod 

QpOD VOD RdS 

In comparison with the maximal disc simulation (run 68) 
from D&S our scaled pattern speed is slightly lower than 
theirs for the duration of the simulation, 0.086 rad t -1 
instead of about 0.1 after the bar is fully formed, and 
0.054 rad t _1 instead of around 0.057 at the end of our sim- 
ulation. Although certainly disc dominated, our mass model 
would be more correctly compared with D&S control run, 
which has a similar disc to halo velocity ratio (their r\ pa- 
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Figure 8. Angular momentum of the disk (middle), halo (bot- 
tom), and total (top). The total angular momentum is well con- 
served and there is a transfer of 10% of the disk angular momen- 
tum to the halo. 
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Figure 9. Determining pattern speed (slope) from luminosity- 
weighted velocity over position. 

rameter). Here, our bar pattern speed is slower initially 
(0.093 rad t -1 instead of around 0.15), but does not slow 
down as quickly or as much as their simulation (final speed 
0.058 rad t~ x instead of around 0.041). Because we are using 
different simulation techniques and mass models, the agree- 
ment here is promising, although it should be noted that our 
pattern speed degradation is consistently lower than D&S. 
This agrees with Athanassoula's (2002b) assertion that bars 
slow down at a faster rate in colder disks. 

We also measured the bar pattern speeds using the ob- 
servational technique described by T&W. Taking observa- 
tions from multiple slits laid parallel to the apparent major 
axis of the tilted disc and finding the luminosity-weighted 
position and average line of sight velocities determines the 
pattern speed: 

n . f .s < W^I > -v , , 

Si„sinm — (2) 

P w < x > -x 

(M&K). This method, although direct, is not simple to 
execute observationally. Since the basis of this relation is 
the continuity equation, common observations using gaseous 
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Figure 10. Four methods of determining bar length: ellipticity of the bar (1 — b/a) - the point where l-b/a is a maximum corresponds 
to the end of the bar; phase angle of the ellipses - where the angle deviates 5dcg; m=2 Fourier amplitude and phase; the bar end is 
measured halfway down the slope on the amplitude plot, and where the phase plot deviates by 5deg. t=150 and t=450 are shown with 
the dashed line indicating the measured bar length for each method. 



emissions are disqualified, as the gas cycle through stars is 
inherently discontinuous. The method is therefore most ac- 
curate for those early-type barred spiral galaxies with the 
least amount of gas and star formation. There are also pre- 
ferred tilt angles to the galaxy to minimize error and diffi- 
culties surrounding determining the line-of-sight velocities. 
It does, however, remove many of the uncertainties inherent 
in the indirect methods. 

Conversely, this technique is relatively simple to per- 
form on simulated data. At various points during the simu- 
lation, snapshots were rotated so the bar was at a 45° angle 
to the x-axis, then tilted about this axis by 30°, similar to 
the observed galaxies listed below. Slits were laid through 
the centre of the nucleus and at ± 1 scale length, parallel to 
the major axis. Fig. 9 plots the luminosity weighted veloc- 
ity versus position, with the assumption of a constant M/L 
ratio. 

This method works very well with the simulated galaxy, 
with speeds agreeing to within around 2 per cent of actual 
for the duration of the simulation when using the known 30° 
inclination of the disc. However, in observational analysis the 
tilt of the galaxy is not known, and these simulations show 
that the outer isophotes of face-on systems are not exactly 
circular. Therefore, determining the tilt of the galaxy by 
assuming this regular shape will lead to errors. Based on the 
distance to NGC 4596, the outermost detectable isophote is 



at 5.8 scale lengths. At t=200 (5.4 Gyr), the axis ratio b/a 
at this radius is 0.72. If assumed circular, the computed 
tilt angle would be 44.2° instead of 30°, which would result 
in a 26 per cent error on the pattern speed as opposed to 
the 2.6 per cent error found with the correct inclination. 
These errors systematically underestimate the pattern speed 
and are reduced as time progresses and the outer isophotes 
become more regular with the disappearance of any spiral 
pattern or time-dependent instabilities. 

The definition of a fast or slow bar depends on the 
length of the bar and the positions of the resonances. The 
former is a point of contention. A&M have shown that any 
one method of determining bar length will not work equally 
well for the duration of the galaxy's lifetime. In this paper 
we use several methods to determine the length of the bar: 
the radius where the phase angle of the elliptical isophotes 
twists, the m=2 component of the Fourier amplitude and the 
m=2 phase component. The latter two are used by D&S, and 
all three are used observationally by Aguerri et al. (2003). 
Originally we had attempted to use Abraham et al.'s (1999) 
description for finding the bar length by choosing the ellip- 
tical isophote with the smallest b/a axis ratio (option (i) 
or L b / a in A&M's comparison of bar length determination 
methods), but we found this method consistently underrep- 
resented the bar in comparison to the other methods. A large 
drop in this measure during the last third of the simulation 
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causes the difference to jump from around 25 per cent to 60 
per cent shorter, although we find similar ellipticity plots as 
A&M for their MD model, where the 1 — b/a measure reaches 
a peak before the end of the bar, then gradually declines. 
The peak occurs at a smaller radius for the later time steps. 
The first three methods and the resulting bar length deter- 
mination from each are plotted in Fig. 10. Fig. 11 shows the 
average bar length using the first three methods as a func- 
tion of time, with error bars corresponding to the maximum 
and minimum lengths calculated. The bar length appears to 
be constant or slightly increasing with time. Also plotted are 
the corotation to bar length ratios, which are all under 1.5, 
indicating the bar is fast throughout the galaxy's evolution. 
D&S quote a final corotation to bar length ratio (Dl : at) of 
1.57±0.27 which they call 'acceptably fast ... only barely so'. 
There is some agreement between our results and theirs, es- 
pecially if we only take the Fourier component methods into 
account, where we find a final ratio of 1.44 ± 0.1. 



4.1 Comparison with observations 

The actual pattern speed of the bar is scaled to the dimen- 
sions of real galaxies and compared with observed pattern 
speeds listed in Table 1. 

The T&W method directly measures pattern speed 
from observables, and is thus more accurate than the other, 
indirect, methods. It has thus far only been applied to a 
handful of galaxies, most of which are listed in Table 1. The 
T&W method depends on continuity, which the life-cycle of 
gas inherently contradicts, so reliable use of this measure- 
ment is restricted to early-type galaxies. 

The scaled pattern speeds found in this simulation are 
initially too slow for all but one galaxy (although two more 
are just barely within the lower bound of the observed er- 
ror bars). However, even though the pattern speed decreases 
with time, the redistribution of the disk particles by the bar 
causes the scale length of the disk to increase with time as 
well. Measured with a double exponential fit to the surface 
density profile, this increase in the disk scale length is fairly 
rapid initially, jumping to 2.3 from 1.0 between t — 50 and 
t = 110, and slowly increasing to 2.8 by the end of the sim- 
ulation. This increase more than offsets the pattern speed 
slow-down, with only one comparison still too low at t — 110 
while two are higher than observed, and the rest are within 
the accepted ranges. By the end of the simulation three are 
again too low, only one is too high, and the rest are ac- 
ceptable (with two at the lower limit of the error bars, and 
two well within the measured range). Therefore, for a con- 
siderable duration, these long-lived bar is rotating at speeds 
comparable to those observed. This reiterates the corotation 
to bar length ration finding above: both the simulated and 
observed bars are "fast" (ratio less than 1.5). 



5 EDGE-ON KINEMATICS 

Some attention has recently been placed on the identification 
of edge-on barred galaxies. K&M, Athanassoula & Bureau 
(1999) (A&B), and B&A have studied the line-of-sight ve- 
locity profile of peanut-shaped, edge-on galaxies and have 
shown that their kinematics are explained by the presence 
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Figure 11. Bar length (triangles) determined by the average of 
three techniques plotted in Fig. 10 (ellipse angle, m=2 Fourier 
amplitude and phase; ellipticity was found not to be a good in- 
dicator in this case). Corotation to bar length ratio is noted with 
squares at the bottom of the plot; the ratio is always below 1.5, 
indicating a fast bar. 



of a bar. The characteristic plot of position versus veloc- 
ity gives a unique figure-eight pattern for barred galaxies, 
since some orbits are depleted in these systems. Our simu- 
lation results have been similarly plotted. After initial bar 
formation, the plot appears similar to the unbarred plots; 
as the bar develops, the central positions darken while the 
density increases, and the areas immediately to the outside 
of this central pole suffer some depopulation, similar to the 
observational plots in K&M and A&B. With bar buckling, 
though, comes increased scatter in the plot, removing any 
forming pattern, and once again leaving a plot which would 
be considered unbarred if viewed observationally. Plots at 
t~50 and t=76 are shown in Fig. 12. 

The original simulation is composed only of a disc and 
halo and was constructed to be in equilibrium with a Toomre 
Q greater than 1. This generates a higher velocity disper- 
sion than would be observed: 0.23 (times maximum velocity) 
tangential and 0.215 radial dispersion at 2.4 scale lengths. 
The radial dispersion increases to 0.258 by t=76. This blurs 
the figure-eight pattern in the K&M plot. Since all the ob- 
servational results used cool gas to determine velocities, our 
results would be much closer to the analytically predicted 
orbits, and the distinctive figure-eight plot. 



6 BULGE MODEL 

To counter part of the possible velocity dispersion problem 
stated above, a model containing a small, compact bulge was 
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Table 1. Comparison of simulated and observed pattern speeds. 



Galaxy 


Ro 


Velocity 


obs. n p 


Sim Q p (t=26) 


Sim fl p * (t=110) 


Sim fl p * (t=500) 




(kpc) 


(km s -1 ) 


(km s -1 kpc -1 ) 


(km s -1 kpc -1 ) 


(km s -1 kpc -1 ) 


(km s" 1 kpc -1 ) 


NGC 4596 (SBa) a 


3.2 b 


120 b 


52 ± 13 


16.4 


27.38 


21.6 


NGC 936 (SB0) c 


3.5 d 


280 d 


60 ± 14 


35 


58.4 


46 


NGC 1023 (SB0) e 


2.9 


270 


87 ±30 


40.7 


68 


53.6 


ESO 139-G009 (SABO)^ 


4.19 


314 


61 ± 17 


33 


55 


44 


IC 874 (SBO)-f 


1.88 s 


187 


41.6 ± 14.3 


43 


73 


57 


NGC 1308 (SBO)-f 


3.6» 


347 


99.4 ±34.8 


43 


71 


56 


NGC 1440 (SB0)f 


1.7 9 


283 


83 ± 10 


74 


124 


98 


NGC 3412 (SBO)-f 


2.29 


205 


57 ± 16 


42 


69 


55 



*Determined with instantaneous scale length (at t=110, i? e =2.3 and at t=500, R e =2.8) 

a Gcrssen ct al. 1999 

b Galaxy properties from Kent 1990 

c Mcrrificld & Kuijkcn 1995 

d Galaxy properties from Kent 1989 

e Debattista et al. 2002 

f Agucrri ct al. 2003 

f Debattista, private communication 
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Figure 12. Line-of-sight velocity distributions from t=50 and t=76. A slight figure-eight pattern is emerging at t=50, but is wiped out 
by the buckling phase at t=76. 



constructed in the lower resolution range (500K particles), 
keeping a similar halo and rotation curve (see Section 2 for 
a more complete description). The evolution of the disc is 
similar to the earlier model in that the bar is initiated a 
little before t=20 (540 Myr); it differs in that the initial 
bar is about 20 per cent smaller, with dumbbell-shaped in- 
ner isophotes when viewed face-on. The bar remains smaller 
than the the previous model, by about the same amount for 
the duration of the simulation. This is in direct contrast with 
the findings of A&M who found their bulge model to have 
a longer bar, although the shape is similar. However, their 
bulge is more massive in comparison to the disk, and more 
extended than ours. Fig. 13 shows t=26 to 250. The initial 
velocity dispersion (0.146 versus 0.23 above for tangential 
and 0.18 versus 0.215 for radial) was lower than our original 
simulation; however, the line-of-sight velocity distribution 
(LOSVD) plots remained similar. The high resolution simu- 
lation showed much more detail and followed trends toward 
the figure-eight shape not seen in the lower resolution plots 



for the same initial conditions. Creating a simulation of the 
bulge mass model with higher resolution, with its lowered 
velocity dispersion, should yield clearer LOSVD plots than 
the original high resolution run, therefore we cannot make 
any conclusions as to whether or not the inflated velocity 
dispersion in the original simulation is responsible for mask- 
ing the figure-eight pattern. 

The pattern speed of the bulge model is higher than 
the original for the duration of the simulation, is reasonably 
steady after t=100 and can be compared with observations, 
as earlier. The result is that the pattern speed scales up to 
an acceptable rate for NGC 1023, ISO 139-G009, and NGC 
1308, is too low for NGC 4596, and too high for NGC 936, 
IC 874, NGC 1440, and NGC 3412 for the whole run. The 
final measured pattern speed for the bulge model (at t=250) 
is the same as that for £=50 in the original disc±halo model. 
Some speed increase is to be expected with the higher cen- 
tral density, as was mentioned in D&S; however, too much 
central density means the disc is no longer unstable to bar 
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Figure 13. Non-contoured logscale face-on views of the bulge model for t=26, 50, 100, and 250. Only the disk particles are shown. 



formation. The bar was even closer to corotation, with a 
Di/as of near one for most of the run, ending with 1.3±§ ;f f 



7 HALO BAR 

A slightly elongated shape was excited in our halo, with an 
axis ratio of around 0.88. If the alignment of this halo bar 
were to coincide with the disc bar, the dynamical friction 
should be reduced. Although the Q p deceleration slows after 
about t=150, suggesting a decrease in torque on the bar, the 
halo bar takes much longer to line up with the disc bar. The 
bar orientations are mostly within 10° of one another after 
t=450, although not locked in at the same relative angle. 
D&S also found an m=2 response in the halo which slowly 
aligns with the disc bar. 



The pattern speed degradation during the later time 
steps may not be solely due to the halo friction. Misiriotis 
& Athanassoula (2000) have shown that the Q p is slower 
in thicker discs. As the bar buckles and the simulation pro- 
gresses, the scale height of the disc increases. This should 
then also slow the bar, even if the halo bar is locked into 
synchronous rotation with the disc bar, and the dynamical 
friction between the disc and halo is at a minimum. 



8 COMPARISON OF EARLY AND LATE 
TYPE BARS 

Elmegreen & Elmegreen (1985) studied the kinematic and 
dynamic properties of different Hubble- types. They found 
distinct photometric differences between early (SBO-SBbc) 
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Figure 14. Surface density (luminosity) for t=0 (dotted), t=50 
(dashed), t=250 (dash-dot), and t=500 (solid). t=50 shows an 
initial flat-topped profile, but the flat section is about half the 
length of the bar. The later times show an evolution to a double 
exponential. 



and late-type galaxies. Early-type galaxy bars appear to end 
at or near corotation and are twice the length of the turnover 
radius of the rotation curve. Their luminosity profile is flat 
compared with the exponential decrease of the rest of the 
disc. Late-type bars are shorter, seem to end near the in- 
ner Lindblad resonance, and are shorter than the turnover 
radius. They have an exponential luminosity profile with a 
different slope than that of the disc. 

The properties of the bar evolves during our simulation. 
By assuming a constant M/L, the density profile of the bar 
can be taken as representative of a luminosity profile (Fig. 
14). Soon after the bar forms, a small section of the bar 
near the core has a fiat density profile, with a near expo- 
nential decrease for the remainder of the bar. Later time 
steps show only exponential profiles. The bar length is near 
the corotation resonance for the duration of the run, and 
the bars are always longer than the turnover radius and 
the inner Lindblad resonance. Fig. 15 contains plots of the 
Lindblad resonances for the bar: Q, Q — k/2, and fi + k/2, 
which are determined based on an azimuthally averaged re- 
distribution of the particles to make the disk axisymmetric. 
The bar length criteria agree with an early-type galaxy, but 
the density profile is late-type. Another point to note is that 
early-type galaxies are also characterized by a large bulge to 
disk ratio, whereas our models, at most, have a very small 
bulge. 

The main results we are comparing in this study (evolu- 
tion of pattern speed, corotation to bar length ratio) may ac- 
tually be similar for all barred galaxies. Kinematical studies 
have shown that bars should be fast even in late-type galax- 
ies (e.g. due to dust lane placement Athanassoula 1992, and 
hydrodynamical simulations Weiner et al. 2001). Therefore, 
our small or no bulge galaxy with fast bar with an exponen- 
tial density profile is fully consistent with a real late-type 
galaxy. 



9 DISCUSSION 
9.1 Resolution 

Along with the 20M particle simulation, a 500K particle 
run with the same initial conditions was conducted for com- 
parison. The pattern speed is very similar, although after 
t=350 the lower resolution model slows a little more than 
quoted above (f2 p = 0.18 vs. 0.20). The bar lengths, based 
only on the m — 2 amplitude and phase change, are in gen- 
eral shorter than in the 20M simulation before t=350, and 
longer thereafter, showing a distinct trend to increasing the 
bar length over time. Most of the values are within the error 
bars of Fig. 11, and there is therefore at least in nominal 
agreement between the two simulations. 

The high resolution simulation showed much more de- 
tail in the bar. For example, isophote twists during the for- 
mation and buckling of the bar were clearly seen, whereas 
they were absent in the lower resolution model. The T&W 
method of finding pattern speed was also more accurate with 
the larger number of particles. The 500K model results were 
within 5-10 per cent of the actual pattern speed as opposed 
to around 2 per cent for the high resolution results. 

By the midpoint of the runs, scale height calculations 
were comparable with the original simulation, but at t=100 
the average scale height for the high resolution simulation 
seems to be higher than the 500K model. This is most prob- 
ably due to the fact that the bar in the high resolution sim- 
ulation has already buckled at this point, whereas the low 
resolution model has not; buckling heats the disc. The scale 
height is of some interest, as it shows the interaction of the 
bar with the disc, and an increased scale height has been 
shown to decrease the pattern speed (Misiriotis & Athanas- 
soula 2000). 

Therefore, it can be concluded that the lower resolu- 
tion is capable of outlining many dynamical properties, but 
higher resolution is needed for some detailed comparisons 
with real galaxies. 



9.2 Code Comparison 

Victor Debattista and Jerry Sellwood generously offered the 
use of their mass model initial conditions to run in our sim- 
ulation code for comparison. Their simulation was run us- 
ing a 3-D Cartesian particle-mesh code, with a grid spacing 
of 0.2i?d, and a total particle count of 600K (102K in the 
disc). Here we run the same particles (run 68 in D&S, max- 
imum disc model) with a tree code (Dubinski 1996). We 
initially had problems with our small softening length, until 
we changed this to correspond with the grid resolution in 
their original simulation (equivalent to a Plummer soften- 
ing of 0.4 scale lengths (1.3 kpc) instead of 0.0125 for the 
halo and 0.005 for the disc), as their particles are set up in 
equilibrium with respect to the mesh. They ran many dif- 
ferent initial conditions during their research, and we found 
the Q—1.5 model (unpublished) to evolve similarly to our 
code. We found a pattern speed of 0.095 rad t _1 at t= 1000 
compared to around 0.07 plotted in their paper (their units), 
and a shorter bar than the original, although when Debat- 
tista and Sellwood (private communication) ran the same 
model in a higher resolution mesh code this result was cor- 
roborated. The equilibrium of the particles with the original 
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Figure 15. Lindblad resonances and the bar pattern speed for t=0, 100, 200, and 450. In all cases the bar is slow; at f.=450 the bar 
length is around a scale length and within the ILR as predicted for late-type barred spirals by Elmegreen and Elmcgrcen (1985) 



mesh is most likely the reason for the disparity. The higher 
resolution leads to stronger forces in the central region of 
the galaxy, which the initial velocity dispersion does not 
support, causing a slight collapse of the disc particles to the 
centre at the beginning of the simulation. This causes the 
rotation curve to rise more steeply, which leads to a shorter 
bar (Combes & Elmegreen 1993). The shorter bar feels less 
friction from the halo, which results in less pattern speed 
deceleration and the higher pattern speeds found here. 

Run 68 in D&S was actually conducted with a Q=0.05 
model. Due to earlier attempts on our models, we were con- 
cerned about the low Q value and the stability of the disc. 
We were unsuccessful in reproducing the output from this 
model, as we were in our own discs with a Q much below 
1.0. Their results may be due to the way the particles were 
set up in relation to their grid, or the grid method may be 
supplying some artificial stability to the model. 

A copy of our low resolution initial conditions were also 
forwarded back to Debattista and Sellwood. Using a new hy- 
brid grid code with softening of 0.0257?^ and a time step of 



0.025 they found bar lengths between 2.5 and 3, well within 
the error bars of Fig. 11. They measured the strength of 
the bar by calculating the m — 2 amplitude of the whole 
disc as a function of time. The resulting strength is very 
similar to ours for the same resolution, however we find a 
greater bar strength for the 20M run between t = 150 and 
t = 450. The increase may be due to the larger number of 
particles available for angular momentum transfer in that 
model. Although we agreed well until t = 150, they found 
a significantly higher pattern speed thereafter, as their sim- 
ulation did not slow down as much as ours. The reason for 
this is unclear at the moment. 

While initial differences were seen, the comparison be- 
tween these codes has shown most of the results to be re- 
producible. That said, direct code comparison is very much 
needed in this field to determine the influence of method 
versus initial conditions. 



© 2002 RAS, MNRAS 000, 1-14 



14 J. O'Neill and J. Dubinski 



10 CONCLUSIONS AND FUTURE WORK 



REFERENCES 



The three N-body galaxy models created here generally con- 
cur with previous research in structure (face-on and edge-on 
shapes of bars) . 

The bars measured here are all fast, with pattern speeds 
in the range of those observed directly in early-type galax- 
ies. Although the pattern speeds of late-type galaxies, which 
would be a more appropriate comparison to these simu- 
lations, cannot be directly measured, kinematical studies 
have agreed that their bars are also fast (e.g. van Albada & 
Sanders 1982; Hunter et al. 1988; Athanassoula 1992; Lind- 
blad et al. 1996; Weiner et al. 2001). 

The corotation to bar axis ratio in the high resolution 
simulation is less than 1.5, which is designated a fast bar. 
This is marginally consistent with the D&S simulations, al- 
though our models did not require a fully maximal disk to 
remain fast. 

A recen t preprint by Valenzuela & Klypin (astro- 
ph/0201028) outlined a simulation which used 780K par 



tides (model B), Hernquist initial conditions, and an adap- 
tive tree-mesh code that has a higher resolution at the centre 
than the less active outer regions. Their spatial resolution 
in the central regions is slightly lower than ours, and they 
use many fewer particles than our high resolution model. 
Many of their results are in agreement with those presented 
here: namely, comparable angular momentum transfer, the 
central density build up created by the bar and its affect on 
the rotation curve, the lack of alteration of the dark matter 
density distribution, and the double exponential disk density 
profile. We disagree in bar length (theirs are much shorter), 
actual pattern speed (they find a more rapid rotation), and 
pattern speed degradation (their pattern speed is constant). 
Our disc to halo mass ratio and Q parameter are different, 
though, which may account for some of the discrepancies. 

By comparing our 20M particle simulation with a 500K 
particle model, it was determined that the basic dynamics of 
the system can be reproduced with lower resolution (pattern 
speed, overall evolution), but if observational techniques are 
to be used on the results, higher resolution is needed. The 
techniques used to varying degrees of success here were the 
T&W method to find pattern speed and the K&M edge-on 
line-of-sight velocity distribution plots to find bar orbits. 

The results presented here agree with much of the previ- 
ous work done on isolated, bar unstable galaxy simulations. 
We have shown resolution is no longer a large source of er- 
ror. We suggest further investigations should concentrate on 
generating early-type barred spiral galaxies with high reso- 
lution low softening codes. Also, more observational data of 
real pattern speeds, especially for late-type galaxies, would 
be a great asset to the field. 
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